Collective synchronization in spatially extended systems of coupled oscillators with 

random frequencies 
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We study collective behavior of locally coupled limit-cycle oscillators with random intrinsic fre- 
quencies, spatially extended over d-dimensional hypercubic lattices. Phase synchronization as well 
as frequency entrainment are explored analytically in the linear (strong-coupling) regime and nu- 
merically in the nonlinear (weak-coupling) regime. Our analysis shows that the oscillator phases 
are always desynchronized up to d = 4, which implies the lower critical dimension df = A for phase 
synchronization. On the other hand, the oscillators behave collectively in frequency (phase veloc- 
ity) even in three dimensions (d = 3), indicating that the lower critical dimension for frequency 
entrainment is df = 2. Nonlinear effects due to periodic nature of limit-cycle oscillators are found 
to become significant in the weak-coupling regime: So-called runaway oscillators destroy the syn- 
chronized (ordered) phase and there emerges a fully random (disordered) phase. Critical behavior 
near the synchronization transition into the fully random phase is unveiled via numerical investiga- 
tion. Collective behavior of globally-coupled oscillators is also examined and compared with that of 
locally coupled oscillators. 

PACS numbers: 05.45.Xt, 89.75.-k, 05.10.-a 



I. INTRODUCTION 

Various systems in nature have been known to ex- 
hibit remarkable phenomena of collective synchroniza- 
tion, which have attracted much attention in the science 
community. To understand collective synchronization be- 
havior, systems of coupled oscillators have been widely 
considered. One of the simplest and typical models for 
those systems was first introduced by Winfree and 
later refined by Kuramoto and others [2> S 13 j consider- 
ing additional ingredients relevant to reality. In existing 
literature systems of fully coupled oscillators have mostly 
been considered due to its analytical tractability and sim- 
plicity. On the other hand, systems of spatially extended 
oscillators with local couplings have received less atten- 
tion even though they are more realistic and frequently 
observed in nature. In some studies on the locally- 
coupled oscillators 0, |E 8, 9, 10], collective synchro- 
nization, in particular, frequency entrainment has been 
investigated. However, they did not provide a clear an- 
swer as yet to whether the ordered (frequency-entrained) 
phase exists in space dimension d = 2, consequently the 
value of the lower critical dimension, and also to the criti- 
cal property near the frequency-entrainment transition in 
higher dimensions. Phase synchronization has also been 
studied; however, there still remain many fundamental 
questions that are not clearly answered, including the 
lower critical dimension for phase synchronization and 
critical scaling properties. 

Collective synchronization in phase or frequency 
(phase velocity) of infinitely many oscillators may emerge 
via competition between ferromagnetic-type coupling 
and scatteredness of random intrinsic frequencies. In the 
strong-coupling limit, oscillators prefer to behave collec- 



tively, overcoming the randomness of intrinsic phase ve- 
locities. In the weak-coupling limit, each oscillator tends 
to be obedient to its intrinsic frequency and makes the 
system disordered in phase and/or frequency. In low 
space dimensions the system with local couplings may 
even not be partially ordered (synchronized) for arbitrar- 
ily strong coupling strength. 

In this paper, we obtain this dimensional threshold 
for appearance of the ordered phase, i.e. the lower crit- 
ical dimension for both phase synchronization and fre- 
quency entrainment. Also unveiled is the nature of the 
synchronization transitions in higher dimensions. In the 
strong-coupling regime, we linearize the equations of mo- 
tion for the oscillators and investigate analytically their 
collective behavior. An analogy to the surface growth 
problem helps us to probe the synchronization order pa- 
rameter. In the weak-coupling regime, we numerically 
integrate the equations of motion and find, via careful 
finitc-size-scaling analysis, that the synchronized phase 
and the phase synchronization transition appear only for 
d > 5 while the frequency entrainment appears for d> 3. 

This paper consists of five sections: Section II in- 
troduces the system of locally-coupled oscillators on d- 
dimensional hypercubic lattices. In Sec. Ill, phase syn- 
chronization is studied by means of the linear theory and 
numerical integration. Nonlinear effects are interpreted 
in terms of runaway oscillators and the analogy to the 
surface landscape is discussed. We also investigate the 
nature of the phase synchronization in five and six dimen- 
sions. This study of phase synchronization is complemen- 
tary to our recent work . Section IV is devoted to the 
study of the frequency entrainment behavior. Discussed 
are nonlinear effects and complete frequency entrainment 
as well as the nature of the frequency entrainment transi- 
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tion. In Sec. V the globally coupled system is examined 
and compared in detail with the locally coupled system. 
Finally, a brief summary is given in Sec. VI. 



II. LOCALLY COUPLED OSCILLATORS 

A general description of coupled oscillator systems may 
be given by 



dt 



— Fi (Xi ) + Gij (X-i , Xj ) , 



(1) 



where X^ is the vector describing the ith oscillator and 
Gij represents the coupling between oscillators i and j. 
When the coupling is not strong (i.e., for small G), one 
can show that only phases (rather than amplitudes) of os- 
cillators are relevant (see, e.g., Ref. 0)- It is well known 
that, as the coupling becomes sufficiently strong, the am- 
plitude variation may not be neglected and one should 
resort to the complex Ginzburg-Landau (CGL) descrip- 
tion; this is often suitable for describing the synchroniza- 
tion phenomena of identical oscillators, which overcome 
the frequency variations originated from amplitude vari- 
ations. 

In this paper, we do not intend to cover the regime 
where the CGL description is valid. We instead focus on 
the collective synchronization transition displayed by a 
large number of limit-cycle oscillators with random in- 
trinsic frequencies, where only oscillator phases are rele- 
vant fluctuating variables. 

We thus consider the set of equations of motion 



dt 



K sm 



{4>i - (t>j ), 



(2) 



which governs the dynamics of N coupled limit-cycle os- 
cillators located at sites of a d-dimensional hypercubic 
lattice. Here 4>i represents the phase of the ith oscillator 
{i = 1,2,- ••,iV), whereas the first term and the sec- 
ond term on the right-hand side represent the intrinsic 
frequency of the ith oscillator and the local interactions 
between the ith oscillator and its nearest neighbors, the 
set of which is denoted by A^, respectively. The intrinsic 
frequency oJi is assumed to be randomly distributed ac- 
cording to the Gaussian distribution function g{ijj) with 
mean which we set = without loss of generality, 
and variance 2a. The coupling is assumed to be fer- 
romagnetic, i.e., X > 0, so that neighboring oscillators 
favor their phase difference minimized. The sine function 
form is the most general representation of the coupling 
in the lowest order and its periodic nature is generic in 
limit-cycle oscillator systems. Higher-order terms are ir- 
relevant in the sense of universality. It is also noteworthy 
that Eq. Q with spatio-temporal random noise (thermal 
noise) rji (t) instead of the quenched noise Ui describes the 
dynamics of the well-known 0(2)-symmetric XY model. 

When the coupling is absent {K = 0), each limit- 
cycle oscillator evolves with its own intrinsic frequency 



uji according to d(j)i/dt = Ui, resulting in that the sys- 
tem becomes trivially desynchronized. For finite cou- 
pling {K > 0), locally ordered (synchronized) regions 
emerge, inside of which oscillators evolve with a coupling- 
modified effective frequency. Here the dispersion of in- 
trinsic frequencies competes against the coupling and this 
competition sets the size of locally ordered regions. When 
the coupling is strong enough to overcome the disper- 
sion of intrinsic frequencies and subsequently to create 
a globally ordered region, the system exhibits collective 
synchronization behavior. 

For the system of coupled limit-cycle oscillators, two 
kinds of collective synchronization may be considered: 
frequency synchronization and phase synchronization. 
The former is often called frequency entrainment or phase 
locking. 



III. PHASE SYNCHRONIZATION 

We first investigate the phase synchronization, which 
may be probed by the conventional phase order parame- 
ter 




N 



(3) 



with (• • •) denoting the average over different realizations 
of intrinsic frequencies. A nonzero value of A then im- 
plies the emergence of phase synchronization. 

Up to date, analytic treatment has been available only 
at the mean-field (MF) level (see Sec. VII). Namely, in 
the case of globally coupled oscillators where each oscil- 
lator is coupled to every other one with equal strength 
K/N, it is known that the system exhibits collective 
synchronization in phase, which is described by A ~ 
{K - K^Y with 13 = 1/2 near the critical coupling 
strength Kc — 2/'Kg{0) Q- It is of interest to note that 
both the phase synchronization and the frequency en- 
trainment emerge simultaneously at the same coupling 
strength in the globally coupled system. 



A. Linear theory 

Systems with locally coupled oscillators have been lit- 
tle investigated, presumably due to the difficulty in ana- 
lytical treatment. The nonlinear nature of the sine cou- 
pling term in Eq. ((SJ is the major obstacle toward ana- 
lytic treatment. Accordingly, we first suppose that, for 
sufficiently strong coupling strength K , the phase differ- 
ence between any nearest neighboring oscillators is small 
enough to allow the expansion of the sine function in the 
linear regime. Taking the appropriate continuum limit 
in space, we obtain Eq. in the form 



d 



(x, t) = w(x) + KV'^cj) + O (VV, (V0)2V20) , (4) 
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where w(x) are uncorrelated random variables, satisfy- 
ing (w(x)) = and (a;(x)cj(x')) = 2aS(x — x'). We note 
that the couphng constant K and the variance 2a may 
be renormahzed through a coarse graining procedure in 
taking the continuum hmit. However, the renormahzing 
factor may be absorbed into the time scale via proper 
rescaling of time. For convenience, we also relax the con- 
straint < (j) < 2'K and extend the range to (—00,00), 
which is permissible in the linear regime. 

Equation with the irrelevant higher-order terms 
neglected, is reminiscent of the well-known Edwards- 
Wilkinson (EW) equation 0|. The EW equation tra- 
ditionally describing surface evolution becomes identi- 
cal to Eq. by interpreting the phase t) as the 
front height of the growing surface. Note, however, that 
the noise lu{x.) is generated not by conventional spatio- 
temporal disorder but by so-called columnar disorder 
which has only spatial dependence. In other words, the 
disorder w(x) is independent of time, which is differ- 
ent from the conventional thermal noise. The coupling 
strength K plays the role of surface tension in the grow- 
ing surface model. 

A key quantity of interest in the context of surface 
growth models is the surface fluctuation width W defined 

by 



(5) 



where L is the linear size of the d-dimensional lattice 
(L"* — N) and ^(t) is the spatial average 



(6) 



Since Eq. I^J is linear, it is convenient to consider the 
Fourier transform 



in terms of which Eq. reads 



d 



(7) 



(8) 



with higher-order terms neglected. The solution of 
Eq. (|SJ) is easily obtained as 



which in turn gives the mean-square width 



1 



d'^k / d'^k {(l){k,t)(f>{k ,t)) 



(27r)2'i 

X (1-6-^'='*) (1-6-^'="*) 

dkk'^-'^il - 26-^*^'* + e-^^'^'^'flO) 



with Qd = 2^-'^'K-'^/'^ /V{d/2) and a denoting the lattice 
constant. In obtaining Eq. (|10|l . we have used the relation 
(u;(k)a;(k')) = 2cr(27r)''J''(k k') and taken 0(k, 0) = 
for the initial condition. 

In the long-time limit [Kt ^ L^), the surface width in 
the stationary state scales, for large L: 



4-d 



~ (CT/47r2i^2)lni, 
~ 2<j/K^, 



d < 4 
d = 4 
d>A. 



(11) 



At any finite values of K , the surface width W thus di- 
verges as L ^ 00 for d < A whereas it remains finite for 
d > A. This indicates that the surface is always rough 
(except dX K = 00) for d < A and always smooth (except 
&i K = Q) for d> A. In the short-time regime {Kt L^), 
in Eq. itinjl becomes 



- {2<j/zK^){Kt)^^-'''>/', d<A 
~ {a/A7r^zK^)\n{Kt), d = 4, 



(12) 



with the dynamic exponent z = 2. For d > A, the mean- 
square width thus saturates to ~ 2a /K"^ exponen- 
tially. 

As expected, the exponents a and /3, which charac- 
terize the width fluctuations according to ~ and 
W ^ \a. the long- and short-time limits, respectively, 
are different from those of the original EW model with 
conventional thermal noise. The difference is attributed 
to the quenched columnar noise which does not have time 
dependence. In fact, simple power counting yields the 
values of a and /3, which are consistent with Eqs. 1)11(1 
and lfn|l . 

We can also derive analytically the exact stationary- 
state probability distribution P [{(/>}] in the linear regime 
described by Eq. (0J. Note that it is usually quite dif- 
ficult to find the exact distribution function for a sys- 
tem with quenched disorder. Equation JSJ assures that 
(/)k should become oj]/^/ {Kli^) for any k in the stationary 
state, where 4>\^ and Wk are used in place of (?!)(k) and Ci;(k), 
respectively, for brevity. Accordingly, one can write the 
stationary-state probability of finding the configuration 
= {0ki, • • • , ^kis,} for given distribution {wk}: 



Averaging over random frequencies {^k}, we find 



(13) 



^ j I?a;kP.,({0k}) 



= / dcjki • • • dwkN exp 
if2 



^k, 
^ Aa 



^kj 



2tt/L 



exp 



exp 



Aa 



i V^ki 



-^jiy^fd. 



(14) 
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The mean-square width VF^ in the stationary state, given 
by Eq. may be derived directly from this distribu- 

tion function. 

It is of much interest to note that this nonequihbrium 
stationary-state probabihty distribution is identical to 
that of the equilibrium Laplacian (tensionless) roughen- 
ing model 13]. The prefactor K"^ /[A:a) in the exponent 
of Eq. H14|l plays the role of the inverse temperature in 
the Laplacian roughening model. One may also relate 
our model to the linear molecular beam epitaxy (MBE) 
growth model |14[. the evolution dynamics of which is 
described by 

^</)(x,t) =7;(x,i)-z.V4<^, (15) 

where 77 (x, t) is the spatio-temporal random noise satisfy- 
ing (?7(x,t)) = and {r]{x,t)r]{x' ,t')) ^ 2DS{x- K')5{t - 
t'). It is straightforward to show that the stationary dis- 
tribution of this MBE model is also identical to Eq. I|14|) 
with the correspondence D a and K'^/2. How- 
ever, the dynamic behavior of the MBE model is char- 
acterized by the dynamic exponent z = 4, in contrast to 
z — 2 vci our model. 

Another important observation is the Gaussian prop- 
erty of the stationary distribution, which provides a link 
between the mean-square width W'^ and the phase order 
parameter A. The practical definition of A in Eq. Q is 
not convenient for analytical treatment, and instead we 
use the formal but equivalent definition: A = (e*'^'^"'^^), 
where denotes the spatial average. The Gaussian prop- 
erty of the stationary probability distribution in Eq. I|14l) 
guarantees (e*-^*^"^') = e"^-^ ('^))/2 for an arbitrary function 
/((/)). Therefore, one can easily obtain 

A = exp [-W'^/2\ (16) 

in the stationary state. It should be noted that this re- 
lation is valid only in the linear (strong-coupling) regime 
described by Eq. Q . 

With this relation, we can express our results for the 
width W in the phase synchronization language: The 
oscillators are always desynchronized (disordered; A = 
0) for d < 4 in the thermodynamic limit {N — > cxj), 
whereas they are always partially synchronized (ordered; 
A ^ 0) for d > 4 at any finite coupling, leading to A = 
exp(-Adcr/iv:2) with Ad = Q.di^'^''^ / (d- '^). Accordingly, 
in the framework of the linear theory, there is neither 
phase synchronization-desynchronization transition nor 
complete phase synchronization (A = 1) at any finite 
coupling K in any dimension d. 

Our linear theory is valid in the strong-coupling 
regime; as the weak coupling regime is approached, the 
original (nonlinear) system described by Eq. Q should 
be more disordered than the prediction of the linear 
theory. This establishes that the full nonlinear system 
should also be desynchronized (A = 0) for d < 4 at any 
finite K ^ which implies that the lower critical dimension 
for the phase synchronization may not be less than four: 



df > A. The nature of the desynchronized phase may 
become different from what is expected from the linear 
theory, especially in the weak-coupling regime. At small 
K, the system becomes far more disordered so that the 
phase difference between nearest neighboring oscillators 
can grow large enough to invalidate the expansion of the 
nonlinear sine function exercised in obtaining the linear 
theory in Eq. For d > 4, it is reasonable to expect a 
phase synchronization transition (roughening transition 
in the surface growth language) at a finite value of K. 
Emergence of the desynchronized phase at nonzero K for 
d> A should be attributed solely to nonlinear effects. Of 
course, one may not rule out the possibility of either the 
full destruction of the synchronized phase at any finite 
K or the absence of the desynchronized phase at all. 

Before investigating the full nonlinear system de- 
scribed by Eq. (0) , we examine the self-consistency of our 
linear theory by considering another standard quantity 
in surface growth models, the height-height correlation 
function 

C(x,t) = ([0(x,t)-0(O,i)]'). (17) 

Similar to the case of W^, the correlation function is 
easily obtained by means of the Fourier transform: 

C(x,t) = ^ /d'^fc J^fi_2e-^'='*+e-2^'='*) 
^ ' ' {2^^ J K^k^ V J 

X [1 - cos(k • x)] . (18) 

For small x (= |x|), the term cos(k • x) may be expanded 
as 1 — k'^x'^/2 + O(x^), which finally yields the following 
stationary behavior: 

C{x) - {2cr/K^)x^L^-'^, d<2 

~ {a/2TTK^)x^\nL, d^2 (19) 
- {2a/K^)x^-'^, d>2. 

In the short-time regime {Kt <C L^), on the other hand, 
we have 

C(x,i) - (2o-/zi^2)^2(^^)(2-d)/^^ 

~ {a/2TTzK^)x'^\n{Kt), d^2 (20) 

with z = 2, and expect exponential saturation for d > 2. 

Note that for d < 2 the correlation C(x, <) diverges 
indefinitely with time in the thermodynamic limit (L — > 
00). With x taken as a lattice unit vector, the corre- 
lation function represents mean-square phase difference 
between nearest neighboring oscillators. For later use, we 
define the mean-square nearest-neighbor phase difference 
(step height) averaged over all lattice directions as 

G'iK,t)^{[V^f). (21) 

For d < 2, the average nearest neighbor phase differ- 
ence G{K) in the stationary state is unbounded for any 
finite K in the thermodynamic limit. Since our linear 
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theory is based on the boundedness of \V(f)\ (in expanding 
the nonhnear sine function and dropping off higher-order 
terms), this imphes that for d < 2 there is no range of 
K where the hnear theory apphes. Therefore the nature 
of the desynchronized phase may be characterized not 
by continuous surface landscape expected from the hn- 
ear theory but possibly by ruptured and splitted surface 
landscape, which will be discussed later. In contrast, for 
d > 2, G{K) remains finite even in the stationary state 
and the linear theory is self-consistent at least for large K 
where G{K) ^ 0(1)- For large G{K), the linear theory 
collapses again, which may give rise to the desynchro- 
nized phase of discontinuous surface character. 
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B. Nonlinear regime 

We now investigate the nonlinear effects which appear 
due to the sine coupling in Eq. |5J). In particular, it would 
be most interesting to probe the possibility of emergence 
of the desynchronized phase at finite coupling strength 
{K 7^ 0) in higher dimensions {d > 5) and the nature 
of the synchronization-desynchronization transition. In 
addition, it would also be of interest to understand the 
nature of the desynchronized phase in lower dimensions 
{d < 4) and possibly a rupturing-type phase transition. 

Unlike in the linearized case, the phase difference |V(/)| 
may not be bounded even in a finite system but diverge 
eventually with finite angular velocity. (Recall that in the 
linear theory the phase difference is always bounded by 
the stationary value of G which is finite for finite system 
size L.) Once the intrinsic frequency difference between 
neighboring oscillators is large enough to win over the 
ferromagnetic coupling, the phase difference may grow 
linearly with time, indefinitely even in a finite system. 

In the regime of weak coupling, those runaway oscilla- 
tors with scattered angular velocities dominate and their 
phases become completely random one another. It is easy 
to show that the phase order parameter defined in Eq. (jS)) 
should decay algebraically as A ^ N^^^^ = L^"^/^ when 
the phases of N oscillators take fully random val- 
ues. On the other hand, in the strong-coupling regime, 
where the linear theory applies, A vanishes exponentially 
as A ~ exp[—{a/4:TT^K'^)L] for d = 3 and algebraically 
as A - i-'y/{8^^K^) for d = 4 [see Eqs. (HTJ and lfTC|) ]. 
The oscillator phases in this regime are desynchronized 
(A = 0) as L ^ 00, but they are correlated. The land- 
scape of these phases exhibits a continuous surface even if 
the characteristic width diverges with the system size. In 
the regime dominated by runaway oscillators, the land- 
scape should be very spiky with diverging width, even 
in a finite system. We call this regime the fully ran- 
dom (desynchronized) phase, while the regime where the 
linear theory applies is dubbed the correlated random 
(desynchronized) phase. For d = 3 and 4, one may expect 
the transition between the fully random and the corre- 
lated random phases. Details of this transition will be 
discussed elsewhere. 




100 



FIG. 1: Phase order parameter A and the mean-square width 
for (a) K = 0.5 and (b) A" = 0.1 in d = 5. 



We check numerically the presence of these runaway 
oscillators in the weak-coupling regime. We integrate 
numerically the full nonlinear equation (|2Jl and measure 
both A in Eq. © and W'^ in Eq. (jTUl). In the hnear 
regime, A and should satisfy the relation in Eq. Hl()|) 
in the stationary state. On the other hand, in the nonlin- 
ear regime with runaway oscillators, Eq. H16|l is no longer 
valid and W is expected to grow linearly with time, with- 
out saturation while the order parameter A should satu- 
rate in a finite system. 

In Fig. n we plot the time dependence of the mean- 
square width and —2 In A at two different coupling 
strengths {K ~ 0.5 and 0.1) in d = 5. Figure^^a) mani- 
fests that for the strong coupling {K ~ 0.5) the relation 
in Eq. (|16|l is well satisfied in the saturated regime (i.e., 
in the stationary state). On the other hand, when the 
coupling is weak {K — 0.1), the breakdown of the rela- 
tion is evident in Fig.^b). In addition, we find that W 
grows linearly with time t in the long-time limit. This 
linear growth starts rather randomly in time, depend- 
ing on the disorder realization (i.e., the distribution of 
random intrinsic frequencies) and also on the initial con- 
dition. Averaging over the disorder realization and over 
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the initial condition, we still obtain the linear growth of 
the width W. Finite-size dependence A ~ N~^^^ also 
supports the presence of runaway oscillators, which will 
be discussed in the next subsection. 



C. Numerical results 

In this subsection, we explore collective phase synchro- 
nization of the coupled oscillators described by Eq. jSJ. 
We integrate numerically Eq. (|2Jl and measure the phase 
order parameter A at various values of the coupling 
strength K and the system size L. For convenience, the 
Gaussian distribution with unit variance {2a = 1) has 
been chosen for the distribution of intrinsic frequencies, 
g{u!) ^ exp(— a;^/4cr), and periodic boundary conditions 
are employed. We begin with the uniform initial con- 
dition {(pi = 0) for given set of {uJi}, chosen randomly 
according to g{uj)- We then use Heun's method [l^ to 
integrate Eq. Q with various values of the discrete time 
increment St. Typically, the integration has been per- 
formed up to TVt = 4 X 10'' time steps with St = 0.05. We 
measure the order parameter A averaged over the data in 
the stationary state, reached after appropriate transient 
time {Kt 3> L'^)- For this purpose, we have discarded 
data typically to the first 2.8 x 10^ time steps. Both St 
and Nt have been varied to check possible systematic de- 
viations; it has been confirmed that so such deviations 
were observed within statistical errors. For the disorder 
average over the distribution g{uj), we have carried out 
one hundred independent runs with randomly chosen re- 
alizations {wi}, over which the averages are taken. 

Figure |21 displays the behavior of the phase order pa- 
rameter A as a function of exp{—K) for d = 2 to 5. For 
d = 2 and 3 [Figs. I2a) and (b)], it is evident that the 
phase order parameter A decreases rapidly as the sys- 
tem size L is increased. It appears to approach zero in 
the thermodynamic limit for any finite K. On the other 
hand, the size dependence of the phase order parameter 
for d = 4 and 5 [see Figs. |2Ic) and (d)] is very different 
from that for c? = 2 or 3. One is thus tempted to conclude 
that A may approach a nonzero value in the thermody- 
namic limit for large K, where the synchronized phase 
emerges. 

We analyze our data in detail by means of finite-size 
scaling and show in Fig. the log-log plots of A versus 
at various values of K. For d = 2 in Fig.|3^a), we ob- 
serve that A ~ as expected, which is a characteristic 
of the fully random phase dominated by runaway oscil- 
lators, up to very large values of K [e.g., e"^ « 0.03]. 
We believe that this fully random phase should extend 
to arbitrarily large values of K, as suggested in previous 
subsections. 

For d = 3 in Fig.|3Jb), this fully random phase seems 
to terminate at a finite value of K. Our data are con- 
sistent with the prediction for the fully random phase, 
A ^ L"^/"^, in the weak-coupling regime {K < Kq) 
and with the prediction for the correlated random phase. 
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FIG. 2: Phase order parameter A plotted as a function of 
exp(— AT), where K is the coupling strength, with the linear 
size L varied for d — (a) 2, (b) 3, (c) 4, and (d) 5. 
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A ~ exp[—{(j/ATr^K'^)L], in the strong-coupling regime 
(K > Kq) where the Hnear theory apphes. The transition 
point between the fully random phase to the correlated 
random one may be determined by the stability analy- 
sis of the linear theory. The average nearest-neighbor 
phase difference G{K) defined by Eq. can be used 
to determine the coupling strength Kq at the transition 
point according to G{Kq) w 0{1). Numerically, we find 
that Kq « j-K, which is equivalent to the analytical 
estimate by setting G'^[Kq) = 1/2. However, in order 
to firmly establish the correlated random phase in the 
thermodynamic limit and explore nature of the transi- 
tion into the fully random phase, one needs much more 
extensive numerical simulations as well as possibly per- 
turbation theory (or stability analysis of the fully random 
phase) in the weak-coupling limit, which are beyond the 
scope of this paper. 

For d = 4, our data in Fig. [2Ic) seem to suggest that 
A remains finite for large K ^ even in the thermodynamic 
limit. However, our analytic argument based on the lin- 
ear theory excludes the possibility of a nonzero value of A 
at any finite value of K in the thermodynamic limit. In 
order to resolve this apparent puzzle, we analyze our data 
carefully by means of finite-size scaling in Fig.|3Jc). Man- 
ifested for K < 0.28 is the fully random phase: A ^ . 
For K > 0.40, A still decreases algebraically with L, as 
shown in the inset of Fig.l^c): A ^ L~'^^-^\ It is pleas- 
ing that our data for K ^ 0.40 agree perfectly with the 
prediction of the linear theory, S{K) — a/Sn'^K'^ (see the 
previous subsection). This result confirms that there is 
no synchronized phase at any finite K for d = 4. It would 
also be interesting to explore the possibility of a phase 
transition near K « Kq = \fo]\ ~ 0.35 between the 
fully random phase and the critical phase described by 
the linear theory; this is currently under investigation. 

For d = 5, it looks evident in Fig. Eld) that there 
exists an ordered (synchronized) phase extended to fi- 
nite values of K. The log-log plot of A versus L^^ in 
Fig. Old) confirms clearly the existence of the synchro- 
nization phase transition. For K < 0.19, we find the 
fully random phase: A ~ L'^/"^. For K > 0.21, on 
the other hand, the inset of Fig. I^d) demonstrates that 
A, first decreasing slightly with L, eventually converges 
to a non-zero value. In fact, for K > 0.24, this sat- 
urated value coincides perfectly with the linear-theory 
value: A = exp[— (T/(127r^iir^)]. Note here that the lin- 
ear theory breaks down for K < Kq — \/cr/9 « 0.24 and 
the transition into the fully random phase apparently oc- 
curs a little later at Kc « 0.20. 

The stable ordered (synchronized) phase begins to 
emerge at d = 5, while the case d = 4 is marginal, appar- 
ently displaying the critical phase; this suggests that the 
lower critical dimension for the phase synchronization is 
— 4. We also note that no complete phase synchro- 
nization (A = 1) is observed at any finite K at least up 
to d = 6 (see Fig.E)). 



FIG. 3: Log- log plot of the phase order parameter A versus 
the inverse size at various values of the coupling strength 
K for d = (a) 2, (b) 3, (c) 4, and (d) 5. Detailed behaviors are 
shown in the insets. 
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D. Phase synchronization transitions for 5d and 6d 

In this subsection we investigate the nature of the 
phase synchronization transition in five {d = 5) and six 
dimensions (d = 6). We numerically estimate the values 
of the critical exponents which characterize the univer- 
sality class of the phase transition. In particular, much 
attention is paid to obtaining the critical exponents /3 
and ly which describe the critical behavior of the order 
parameter and the correlation length, respectively: 

A^{K-K,f and ^-^\K-K,\-^, (22) 

where is the critical coupling strength at the transi- 
tion. 

In a finite system, we assume the finite-size scaling 
relation 

A = L-/5/-/[(if -ifjLi/''], (23) 

where the scaling function behaves f{x) ^ a.s x ^ 
+00 and f{x) ^ constant as a; —> 0. At criticality {K = 
Kc), it reduces to 

A{K,,L)^L-^/r (24) 

To estimate efficiently the exponent fi/v and the transi- 
tion point Kc, we introduce the size-dependent effective 
exponent 

/3 _ ln[A(£0/A(L)] 
i/(L) ln(L'/L) ' ^ ^ 

where we take L' = L + I here. In the thermodynamic 
limit the value of the effective exponent is expected to 
approach zero for the ordered phase {K > Kc), P/v at 
the transition {K = Kc), and d/2 for the fully random 
phase (K < Kc), respectively. 

In Fig. 21 we plot the effective exponents computed at 
various values of K versus L^^ for d = 5. The data for 
K ^ 0.19 are observed to converge to the weak-coupling 
value 5/2, while those for K ^ 0.21 converge to zero 
within statistical errors. Only the data aX K = 0.20 ap- 
pear to converge to a nontrivial value. We thus estimate 
the critical couphng strength Kc — 0.200(5) and the ex- 
ponent ratio P/iy — 1.6(3). 

To check the finite-size scaling relation directly, we plot 
ALl^/" versus {K/Kc-l)L^/'' in Fig.Hand find that the 
data for various values of L and K are collapsed best to 
the curve with the choice Kc = 0.200(5), /S/v = 1.4(3) 
and V = 0.45(10), which results in /3 = 0.63(20). As 
expected, the resulting scaling function f{x) converges 
to a constant for small x, and diverges as x^ for large x 
(see Fig. EJ. Our results for d = 5 arc thus summarized: 

= 1.5(3), J/ = 0.45(10), i^c = 0.200(5). (26) 

We note that there apparently exist substantial devia- 
tions from the mean-field (MF) values, (3/v = \ and 
V = 1/2, although the latter may not be totally excluded. 
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FIG. 4: Effective exponent P/u{L) versus L ^ for d = 5 at 
various values of K. 
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FIG. 5: Data collapse of AL'^Z" against {K/ Kc-1)L^/'' in the 
log-log scale for various values of the system size and coupling 
strength. The best collapse is achieved with P/v = 1-4(3) and 
V = 0.45(10). The straight line has the slope 0.63, giving an 
estimation of [5. 



In view of the argument for the MF nature, these appar- 
ent deviations are rather unexpected and their origin is 
unclear at this stage. 

Similarly, Fig. |B1 displays the plot of the effective ex- 
ponent P/v{L) versus for d = 6, which leads to the 
estimation: 

/?/i/ = 1.0(3), = 0.45(10), ifc = 0.156(2). (27) 

These exponent values appear consistent with the MF 
values, but due to rather large statistical errors it may 
not be conclusive that the d = 6 system exhibits MF-type 
critical behavior. Further investigations are requisite for 
determining unambiguously the upper critical dimension 
for phase synchronization. 
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IV. FREQUENCY ENTRAINMENT 

In this section, we explore the frequency entrainnicnt, 
which is another kind of synchronization behavior ap- 
pearing in coupled oscillator systems. At zero coupling 
{K = 0), the phase of each oscillator evolves with its 
intrinsic frequency uf. d(f>i/dt — tOi. Hence, the fre- 
quency (phase velocity) distribution of the oscillators re- 
mains unchanged from the initial random distribution 
g{uj)- With finite local ferromagnetic coupling {K > 0), 
the oscillators may tend to form locally ordered regions 
where they evolve with an identical frequency. The com- 
petition between the ferromagnetic coupling and the dis- 
persion of intrinsic frequencies sets the size of the locally 
ordered regions. For sufficiently strong coupling, the lo- 
cally ordered regions may expand and merge together 
into a globally ordered region and the number of oscil- 
lators with an identical frequency become macroscopic. 
In this case, frequencies are said to be entrained macro- 
scopically. 

In characterizing the collective behavior of oscillators 
in frequency, two different frequency order parameters 
are introduced. In most previous (numerical) studies 0, 
1^1 1^ 1^ I a frequency order parameter has been defined 
to be 



r= lim {Ns/N), 



N 



(28) 



where Ng and N are the number of oscillators in the 
largest cluster having an identical frequency and the to- 
tal number of oscillators, respectively. Nonzero values of 
r implies the emergence of a macroscopic cluster of oscil- 
lators with an identical frequency. With this definition, 
the frequency-entrainment transition may be viewed as a 
percolation-type phase transition with a continuous fluc- 
tuating variable (frequency). In order to measure r, we 
need geometric information on how the oscillators with 
an identical frequency are placed on the d-dimensional 
hypercubic lattice. 



In the globally connected oscillator system, there is 
generically no geometric information on the placement 
of oscillators, each being connected to every other os- 
cillator. Therefore, in this case, Ng represents simply 
the maximum number of oscillators with an identical fre- 
quency. For later use, we define another frequency order 
parameter without geometrical information as 



= lim (N,JN) 

TV— >oo 



(29) 



where Nn is the maximum number of oscillators with an 
identical frequency. In the system of globally connected 
oscillators, both definitions are equivalent: r = Q. Near 
the frequency-entrainment transition, it is known that 
Q - {K-Kcf with /3 = 1/2 and = 2/7r5(0) % No- 
tice that for globally coupled oscillators both the phase 
synchronization and the frequency entrainment transi- 
tion occur at the same coupling strength and their order 
parameter exponents share the same value. 

For the locally connected oscillators, on the other 
hand, these two definitions are not identical. The order 
parameter Q is always larger than r, because the former 
counts additional oscillators with an identical frequency 
belonging to different clusters. We presume that the or- 
der parameter Q should be more suitable for describ- 
ing the frequency-entrainment transition as an order- 
disorder transition. In general, the percolation-type tran- 
sition characterized by r may or may not occur simulta- 
neously at the same coupling strength with the order- 
disordcr-type transition characterized by Q. However, in 
our model with a continuous fluctuating variable, it may 
be reasonable to assume that these two order parameters 
behave similarly, at least qualitatively near the transition 
from the ordered side. 

Here we measure the order parameter Q to probe the 
frequency-entrainment transition for simplicity and con- 
venience. Accordingly, it is not necessary to retain geo- 
metrical information during integration of Eq. Q . 



A. Linear theory 

Similarly to the case of phase synchronization, we be- 
gin with the linearized equation of motion in Eq. Q and 
measure the fluctuation width of the growing velocity 
(rather than the height) of the surface, which would pro- 
vide key information on the dispersion of the phase ve- 
locity (frequency). The mean-square fluctuation width 
for the phase velocity v(x, t) [= (/)(x, t)] is defined to be 



(30) 



where v{t) is the spatial average. 

Taking the time derivative of the Fourier-space solution 
in Eq. (|5J), we find 



i;(k,t) = w(k)e" 



iffc20(k,O)e 



-Kk^t 



(31) 
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One can easily see that, at any finite K, all Fourier com- 
ponents of the phase velocity except the k = mode 
vanish in the long-time limit it cx)), which indicates 
that the phase velocity becomes uniform in space and V 
approaches zero. Without couphng (i^ = 0), the velocity 
distribution should be identical to the initial frequency 
distribution u;(k) so V'^ is equal to 2a at all times. From 
the above equation, it is straightforward to show the fre- 
quency width, normalized to the fully random value: 

^^nj *. (32) 

In the short-time regime {Kt ^ L^), it decays alge- 
braically for any nonzero K in any space dimension d, 

f ^^-'^/^ (33) 

and eventually vanishes in the regime Kt 3> i^, as ex- 
pected. In the language of surface growth, this corre- 
sponds to the completely flat phase. 

The uniform distribution of the phase velocity {V — 0) 
implies complete frequency entrainment. Therefore, our 
linear theory predicts that, in the stationary state, the 
frequencies are completely entrained (with the frequency 
order parameter r = Q = 1) for any nonzero K in any 
space dimension. Only is there a trivial first-order tran- 
sition a.t K — from the fully random phase {r — Q — Q) 
to the completely entrained phase. Of course, any pre- 
diction for d < 2 should be untrustworthy because there 
nonlinear effects dominate in the whole range of K, inval- 
idating the linear theory. On the other hand, for d > 2, 
the linear theory may survive to establish the completely 
entrained phase for large K where the average nearest- 
neighbor phase difference becomes G{K) ^ C(l). 

Finally, we note that, unlike the phase synchroniza- 
tion problem, there is no explicit and quantitative rela- 
tion between the frequency fiuctuation width V and the 
frequency order parameter r or Q. 

B. Nonlinear effects 

The results from the linear theory can be understood 
rather in a simple manner. Consider the linearized equa- 
tion of motion in Eq. without higher-order terms. 
Taking the time derivative of this equation, the disorder 
term a;(x) drop out and yields a simple noise-free equa- 
tion 

^z;(x,i) =i^V2«(x,t). (34) 

This equation can be viewed as the standard diffusion 
equation governing heat conduction with the diffusion 
constant K and the phase velocity (frequency) field iden- 
tified as the temperature field. Since the diffusion con- 
stant K is positive, one can easily expect that any lo- 
cal temperature gradient should disappear, giving rise 



to a uniform distribution of the temperature field in the 
long-time limit. In terms of the oscillator language, the 
phase velocity of all oscillators become identical, lead- 
ing to a delta- function-like distribution (i.e., completely 
entrained phase in frequency). 

Now we come back to the original equation in Eq. ^ to 
accommodate nonlinear effects. After taking the appro- 
priate continuum limit in space and dropping off higher- 
order terms (but not expanding the sine function to re- 
tain the nonlinear effects at least in lower orders), we 
get 

^<^(x, t) = a;(x) + 2if ^ sin Q (e^ • V)^ ^ , (35) 

where denotes the unit lattice vector in the [i direction 
and the summation is over all d different directions. The 
linearized equation in Eq. Q is restored by expanding 
the sine function and keeping the lowest order. With the 
rotational symmetry (if not spontaneously broken), this 
equation can be written approximately in a simpler form 

^0(x,i) =^(x) + 2ifsinQv20^ , (36) 

which, upon taking the time derivative, gives 

^i;(x, t)=K cos Q V20^ V^v. (37) 

This equation can also be viewed as the diffusion (heat 
conduction) equation, with the effective diffusion constant 
K cos{^\/^ (j)) varying in space and time. More impor- 
tantly, it may become negative depending on the value 
of V^0. With negative diffusion constant, a local ther- 
mal gradient does not diminish but increase and the sys- 
tem can become highly inhomogeneous. As the effective 
diffusion constant changes its sign frequently in time and 
also in space, the system may reach, through competition 
between diffusion and localization (negative diffusion), a 
stationary state with a nonuniform temperature distri- 
bution. In the oscillator language, the phase velocity can 
take a broad distribution in addition to or in the absence 
of delta-function-hke peaks. 

The frequency fluctuation width V representing the 
broadness of the frequency distribution is expected to be- 
come nonzero in the strong-coupling regime where non- 
linear effects may become dominant. We integrate the 
full nonlinear discrete equation in Eq. numerically 
and directly measure V^. Figure displays the normal- 
ized mean-square width of the phase velocity V'^ /2a ver- 
sus exp{~K) with 2ct = 1 and L = 25600, 128, 64, 16, and 
8 for d = 1,2,3,4, and 5, respectively. At zero coupling, 
V^/2a should become unity, while in the completely en- 
trained phase it should vanish. We do not find any ap- 
preciable finite-size effects in the whole range of K in all 
dimensions except for very small values of V'^ /2a, which 
will be discussed later. 
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FIG. 7: The mean-square width of the phase velocity, 
which is divided by the variance 2cr (= 1), is displayed as a 
function of exp{—K) in various dimensional systems. The 
detailed behavior near = is shown in the inset. 

Observed are various interesting features: First, we 
find that the fully random phase exists only at zero cou- 
pling in all dimensions. The normalized width V'^/2a 
becomes smaller than unity as soon as the coupling is 
turned on. This result is in sharp contrast to the case 
of phase synchronization, where the fully random phase 
extends over the whole range of K in one and two dimen- 
sions and also forms a sizable region in the weak-coupling 
regime in higher dimensions. Of course, these two results 
are not contradictory. The oscillator phases cannot be 
correlated without correlated frequencies; however, the 
reverse is not necessarily true. With frequencies already 
correlated, the phases can still be fully random if the 
coupling between oscillators are weak. It should not be 
mistaken that V'^/2a < 1 does not guarantee any macro- 
scopic frequency entrainment (i.e., any nonzero value of 
the frequency order parameter Q). In the case of phase 
synchronization, a similar situation has been found for 
d = 3 and d = 4, where phases are correlated but the 
phase order parameter A remains zero. Therefore, it 
seems that the normalized width does not provide proper 
information on the frequency entrainment-detrainment 
transition. 

Second, the normalized width gradually decreases and 
approaches zero as K is increased in all dimensions. The 
one-dimensional case is special in that the normalized 
width appears finite at any finite indicating no com- 
pletely entrained phase, and approaches the K — oo 
point with a nonzero (very large) slope as a function of 
exp{—K) (see the inset of Fig. |7||. On the other hand, 
all other cases d > 2 appear to possess regions of van- 
ishing width {V — 0) in the strong-coupling regime, 
where V'^ /2a decreases exponentially with an infinites- 
imally small slope at finite values of K. This implies 
that the completely entrained phase may be present at 
finite coupling strength for ci > 2, where the linear theory 
applies. 

However, the stability analysis of the completely en- 




FIG. 8: Behavior of /2a of the phase velocity in various 
dimensional systems, depending on the coupling strength K. 
For comparison, the data for the globally coupled system are 
also shown, represented by gl. 



trained phase should be done with great care. The glob- 
ally coupled system, where each oscillator is coupled with 
every other oscillators with equal strength K/N, pro- 
vides analytic results, with which our results may be 
compared. It is well known that there is no complete fre- 
quency entrainment at finite K for the globally coupled 
system, given with the intrinsic frequency distribution 
g{Lu) having no cutoff at finite frequency u [our choice 
g{uj) ~ exp(— a;^/4(T) provides an example; see the next 
section]. Figure |S1 displays the behavior of V"^ /2a versus 
K = Kz for d = 3, where z (= 2d) is the coordination 
number. For comparison, the data for the globally cou- 
pled system of the same size {N — L"*), with K = K are 
also shown. 

As the system size L is increased, the value of K — 
beyond which the mean-square width for the phase veloc- 
ity vanishes {V'^/2a — 0) tends to shift to larger values, 
albeit slowly, both in the globally coupled and the lo- 
cally coupled systems for d = 3. The value of Kj^ should 
diverge to infinity as L — > oo in the globally coupled sys- 
tem. We do not find much difference between these two 
systems in the value of and its finite-size behavior. 
Actually, the value of is a little bit larger (corre- 
sponding to the narrower completely entrained phase) in 
the d = 3 case than in the globally coupled one. Based on 
this observation, we suggest that for a locally coupled 
system also grows arbitrarily large in the thermodynamic 
limit and there may not exist complete frequency entrain- 
ment at finite K in any dimensions. Of course, if the 
distribution 17(0;) has finite cutoffs (like the semi-circle 
distribution), the locally coupled system might have the 
completely entrained phase at finite values of K, simi- 
larly to the globally coupled one. 
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C. Frequency order parameter 

We now measure the frequency order parameter Q de- 
fined in Eq. (|29|l . whicii is the maximum fraction of the 
oscillators having an identical frequency. 

Each oscillator starts to evolve with its intrinsic fre- 
quency given by the initial distribution g{uj), but the 
coupling between near neighboring oscillators will mod- 
ify its frequency continuously with time until the the sta- 
tionary state is reached. We measure the mean frequency 
of the ith oscillator after some transient time ts'. 



LOi = lim 

t — 'OC 



t-t. 



(38) 



We integrate Eq. (0) up to A^f = 4 x 10'' time steps 
with time increment 5t — 0.05 and measure uji at the 
maximum time i„i — Ntdt with the transient time ts — 
2.8 X W^St. Since the frequency resolution is limited by 
5uj = TT^tm — ts)"^ in numerical integration, the frequen- 
cies of different oscillators are regarded as identical if the 
frequency difference does not exceed 5uj ^ 5.2 x 10~^. 
With this resolution, we can draw a histogram /i((D) of the 
number of oscillators with an identical frequency Co nor- 
malized by the total number of oscillators, i.e. h{u)) — 
1. The order parameter Q is then obtained from the 
maximum (peak) value of h{u)). Precisely speaking with 
finite resolution 5lj, Q \s given by the maximum value of 
h{u)) minus g{Q)5uj = 5Loj\/ Ancr w 2.1 x 10~^, which is 
the maximum value at i^T = 0. In practice, we measure 
this quantity, which should vanish in the detrained phase 
even with finite resolution. 

Figure 1^1 displays the behavior of the frequency order 
parameter Q as a function of exp(— i^T) for d = 1 to 5. For 
d = 1 shown in Fig. [Hfa) , it is evident that the frequency 
order parameter Q decays rapidly with the system size L, 
seemingly approaching zero in the thermodynamic limit 
for any finite K. This implies that there is no frequency 
entrainment at all in one dimension. For d — 2, Fig.l^Jb) 
shows that Q decreases slowly and the entrained phase 
continues to shrink with L. For d = 3,4, and 5, it is 
observed in Figs. Etc), (d), and (e) that Q appears to 
saturate to a nonzero value in the strong-coupling regime, 
which suggests that there exists a frequency entrainment- 
detrainment transition for d > 3. Moreover, it appears 
that the fully entrained phase {Q — 1) begins to show 
up at finite values of K. However, as discussed in the 
previous subsection, the careful analysis in comparison 
with the globally coupled system suggests that complete 
entrainment occurs only at K = oo. 

We analyze our data more carefully to locate the tran- 
sition point Kc separating the entrained phase from the 
detrained one for d — 2 and 3. We define the half-order 
value of the coupling strength JChaif at which the fre- 
quency order parameter becomes one half: Q (if half) = 
1/2 and investigate its finite-size behavior. In Fig. llOr a). 
the half-order value i^haif for d = 2 is plotted ver- 
sus in the semi-log scale, which displays logarith- 
mic divergence ivThaif — alnL with a « 0.35(15). This 
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FIG. 9: Frequency order parameter Q plotted as a function 
of exp{—K), with the linear size L varied for d = (a) 1, (b) 2, 
(c) 3, (d) 4, and (e) 5. 
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FIG. 10: Half-order coupling strength isThaif for frequency en- 
trainment in (a) two (d = 2) and (b) three dimensions {d — 3), 
plotted against the inverse linear size in the semi-log scale 
and in the log-log scale, respectively. For d = 2, JsThaif diverges 
logarithmically with L: Khsdi — alnL with a « 0.35(15). 
For d = 3, T^haif seems to converge exponentially to a value 
around 0.73. The inset demonstrates that the successive slope 
cl between neighboring data approaches zero as L ^ oo. 



confirms that there is no entrained phase presented in 
two dimensions. It is of interest to note the behavior 

1/2 

^haif ~ (In L) in the case of phase synchronization for 
d = 4 (see Sec. III). 

In three dimensions, the situation is clearly different 
from that in two dimensions. Figure Efb) displays the 
log-log plot of iChaif versus As L oo, i^haif is 

observed to converge to a finite value around 0.73. The 
analysis on the successive slope between data points 
for subsequent sizes confirm this convergence. There- 
fore, we conclude that there is a frequency entrainment- 
detrainment transition at K = ^ 0.73 in three di- 
mensions. Similar conclusions can be drawn in higher 
dimensions as well. 

For more quantitative analysis on the nature of the 
detrained phase and the entrainment transition, we show 
in Fig. ^] the log-log plot of Q versus at various 
values of K . Similarly to the phase order parameter A, 
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FIG. 11: Frequency order parameter Q plotted as a function 
of the linear size L for various values of the coupling strength 
K and dimension d = (a) 1, (b) 2, (c) 3, (d) 4, and (e) 5. 
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we expect Q ^ N ^ L '^1'^ for the fully random 
phase. 

Figure lTlT a) shows that in one dimension the frequency 
order parameter decays algebraically as Q ^ i where 
c varies very slowly from the value around 0.6 at K = 3.2 
to the value 0.5 at if = 0. We presume that the appar- 
ent variation of the scaling exponent c with K reflects 
small-size effects and that the one-dimensional system of 
sufBciently large size exhibits the fully random phase in 
frequency for any finite if, with the behavior Q ^ L^^/"^. 

In two dimensions, we also face with a similar situa- 
tion, except that c varies from the value around 0.7 at 
K = 1.0 to 1.0 at if = 0.1. We do not rule out the pos- 
sibility to have a correlated random (detrained) phase, 
similar to the three- and four-dimensional phase synchro- 
nization problems. However, both statistical and system- 
atic (small-size) errors hinder us to clarify this point with 
present data. 

In higher dimensions, there is a frequency entrainment 
transition and we estimate the critical coupling strength 
Kc ~ 0.71(2), 0.32(1), and 0.21(1) for d = 3, 4, and 5, 
respectively. In the entrained phase for K > Kc, the or- 
der parameter Q saturates to a nonzero value as L ^ oo, 
while in the detrained phase [K < Kc), it decays to 
zero. We examine nature of the detrained phase and find 
that it has a characteristic of the fully random phase in 
five dimensions: Q ~ L~^-^. We also investigate the 
critical decay, to find Q - L-l^/" with [3/v = 0.35(15). 
For d = 4, we observe that Q ~ L~^-° in the weak- 
coupHng regime (at K « 0.15), but it decays slower as K 
is increased. However, the variation of the exponent is 
rather small (between 1.7 and 2.0), which suggests that 
the entire detrained phase is also the fully random phase 
in four dimensions. The analysis at the critical point 
reveals that (3/v = 0.2(2). The three-dimensional situa- 
tion is marginal and the present data do not provide even 
suggestive information about the nature of the detrained 
phase. 

It is remarkable that the phase synchronization and 
the frequency entrainment transition apparently occur 
simultaneously (within error bars) for d = 5, just like 
in the globally coupled system. It may be conjectured 
that these two synchronization transitions always oc- 
cur simultaneously for all d > b. However, the critical 
exponents for phase synchronization and for frequency 
entrainment are clearly distinct: {(3/i>)p = 1.5(3) and 
{f3/v)F — 0.35(15), respectively. Note that the two criti- 
cal scalings are identical in the globally coupled system. 

Finally, we study the normalized histogram h{u)) of the 
number of oscillators with an identical frequency uj in the 
stationary state. Figure [T^ exhibits the semi-log plot of 
the histogram h[u)) at various values of K for d — 2 and 
3 with size L — 256 and 64, respectively. At if = 0, 
the histogram should be identical to the initial Gaussian 
distribution g{io). As if is increased, one can see easily 
that the distribution becomes narrower, with its peak be- 
coming sharper. At the same time, its tail part becomes 
thinner and decays faster. This manifests that the os- 
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FIG. 12: Semi-log plot of the normalized histogram of the 
stationary-state frequency of oscillators for d = 2 and 3 with 
L — 256 and 64, respectively. The thin curved line represents 
the Gaussian distribution g{uj). 



cillators having high intrinsic frequencies tend to adjust 
themselves to co-evolve with other oscillators having low 
intrinsic frequencies, through couplings between them. 
This distribution is completely distinct from that of the 
globally coupled system (see Fig. ^| in the next section). 

In summary, our numerical analysis shows that the fre- 
quency entrainment transition emerges for d > 3, which 
suggests that the lower critical dimension for frequency 
entrainment: df = 2. For d > 5, the frequency and the 
phase synchronization transitions occur simultaneously 
but with different scaling behaviors. 



V. GLOBALLY COUPLED OSCILLATORS 

For comparison with locally coupled oscillators, we dis- 
cuss some analytic and numerical results for globally cou- 
pled system, where each oscillator is coupled with ev- 
ery other with equal coupling strength K/N. The set of 
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FIG. 13: Phase order parameter A in the globally coupled 
oscillator system shown as a function of exp{—K) for various 
values of N. 
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FIG. 14: Phase order parameter A plotted against A'^"^ in 
the system of globally coupled oscillators, for various values 
of the coupling strength K. 



equations of motion reads 



K 



N 



LUi > sv[v{(hi — (hi). 

dt 



(39) 



i=i 



This system has been much studied j| , partly due to its 
analytical tractability. Similar to the system of locally- 
coupled oscillators, its collective behavior in phase is de- 
scribed by the order parameter A defined in Eq. Q . For 
convenience, we define the complex order parameter ac- 
cording to Ae*^ = Y.j e**' - With the help of A and 
9, Eq. H39|l reduces to N identical decoupled equations 



dt 



KAsm( 



(40) 



where A and 9 are to be determined by imposing self- 
consistency. 

The stationary solution with constant 9 is then ob- 
tained for the symmetric distribution of intrinsic frequen- 
cies g{uj) = g{—oj). A straightforward calculation leads 
to the self-consistency equation 



A = aKA - c{KAf + 0{KAf 



(41) 



with a = (7r/2W0) and c = -(7r/16)g"(0), which yields 
Kc = 2/'Kg{0) [3. When the distribution g{uj) is concave 
at w = 0, i.e., g (0) < 0, the nontrivial solution (A 7^ 0) 
appears via a pitchfork (supercritical) bifurcation as K 
is raised beyond Kc- Near Kc, the nontrivial solution 
behaves as 



A oc (i^ - Kc) 







(42) 



where /3 — 1/2 is the mean- field value of the order pa- 
rameter exponent. 

We have also integrated numerically Eq. (|39|l and show 
in Fig. ^1 the obtained behavior of the phase order pa- 
rameter A depending on exp(—K) for various values of 



N. As expected, there is a phase synchronization tran- 
sition at K = Kc « 1.6 (or e~^^ « 0.20), which is con- 
sistent with the analytic value Kc = 2/'Kg{Q) = -^/S/tt w 
1.596. 

For more quantitative finite-size analysis, we show in 
Fig.^^the log- log plot of A versus A^~-^. For K < Kc, we 
find that A ~ N~^/'^, which implies that the desynchro- 
nizcd phase is the fully random phase. Near K — Kc, we 
assume the finite-size scaling relation 



A = N-PI^'F 



{K - Kc)N 



1/9 



(43) 



where the critical exponent D describes the divergence 
of the correlation volume (the number of correlated os- 
cillators) (^v K Kc- In d dimensions, we expect 
~ ^'^ with the correlation length ^, which leads to 
the relation v = vd. For the globally coupled system, 
neither the length scale nor the space dimension have 
proper meaning, so only may be properly defined. 
The scaling function behaves F[x) ~ x'' as a; ^ 00 and 
F{x) ^ const, as x ^ 0. 

Ai K = Kc, we have A ^ N^^/'^ , which is shown in 
Fig. USl The best fit is obtained with fi/v = 0.21(2), 
which, together with the exact value /3 = 1/2, leads to 
the estimate v — 2.4(2). Summarizing our results on the 
phase synchronization in the globally coupled system, we 
write 



l3=\/2 and z7 = 2.4(2). 



(44) 



Note that our present estimate is somewhat higher than 
the existing one 2.0(2) fl^. We have also checked the 
finite-size scaling relation directly by plotting AN^^'^ ver- 
sus (K/Kc — 1)A^^/'^ and found the consistent value of 
D — 2.4(2). This may provide a hint on the upper crit- 
ical dimension d^ for phase synchronization. Following 
Ref. IT^ one may assume that the relation D = vmf du 
holds in the coupled synchronization problem. Then, 
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FIG. 15: Critical decay of the phase order parameter. The 
straight line represents the best linear fit of slope 0.21. 
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FIG. 16: Semi-log plot of the normalized histogram of the 
stationary-state frequencies of oscillators in the entrained 
phase {K = 1.7). The thin curved line represents the Gaus- 
sian distribution G{ijj) 



with the usual mean field value umf = 1/2, we come 
up with K, b. If this turns out to be correct, our re- 
sult for the locally coupled system in five dimensions may 
reflect the log-type corrections to the mean-field values. 
At present, however, it is too early to conclude on the 
value of the upper critical dimension. 

We now study the frequency entrainment behavior for 
the globally coupled oscillator system. Similarly to the 
case of the locally-coupled system, we investigate the nor- 
malized histogram h{u)) of the number of oscillators with 
an identical frequency Co in the stationary state. Fig- 
ure E| shows the histogram obtained numerically in the 
frequency entrained phase {K > Kc)^ which exhibits a 
delta-function peak at the entrained frequency {uj — 0) 
in the background Gaussian-type distribution. Looking 
carefully at the background distribution, we find two 
symmetric humps near w = and its Gaussian shape for 
large \u)\ is almost identical to the initial frequency dis- 



tribution g{uj) 0,0. On the other hand, in the detrained 
phase {K < Kc), there are neither peak nor humps and 
the histogram h{u)) is identical to g{uj). This implies 
that the detrained phase for the globally coupled system 
is fully random. (Of course, in a finite-size system, the 
frequency order parameter is nonzero at all K, so one can 
see the peak-hump structure even for K < Kc- However, 
in the thermodynamic limit, the peak- hump structure 
should disappear for K < Kc-) 

As K increases beyond Kc, the oscillators with very 
low intrinsic frequencies {uji w 0) are entrained first, so 
the histogram shows two humps and dips near the en- 
trained frequency {u — 0). This distribution is clearly 
distinct from that of the locally coupled system, where 
the oscillators with high intrinsic frequencies are affected 
first due to local coupling with nearby oscillators with 
low intrinsic frequencies (see Fig. ^] and Sec. IV). In 
the globally coupled system, all oscillators are coupled 
with the same strength (no local environment involved), 
and accordingly oscillators with similar frequencies are 
entrained first. 

The analytic form of the normalized histogram in the 
thermodynamic limit is given by 0, |^ 



(45) 

where the frequency order parameter Q is determined by 
the normalization condition. A concise expression for Q 
then reads 



KA 



g{Lu)duj. 



(46) 



KA 



For A = 0, one can easily see that /i([D) — g{Cj). On 
the other hand, for A > 0, a delta- function peak and 
two humps appear in Eq. (|45|l . as expected. For small A 
(near criticality) , the frequency order parameter behaves 
as 

Q ~ 2KAg{0) a: {K ~ Kef with /3 = 1/2. (47) 

Note that both phase synchronization and frequency en- 
trainment share the same critical point and the same 
scaling behavior near criticality. 

We also study the finite-size scaling behavior of Q at 
criticality and show the log-log plot of Q versus N^^ in 
Fig. ini It is observed that Q - iV^'^/p with (3/D = 
0.21(2), which is also consistent with the value for the 
phase order parameter A. In the detrained phase, we 
find Q ~ N~^^'^, as expected. 

As the humps disappear as K one may de- 

fine the frequency distance D between the delta peak 
position (a) — 0) and the hump position, as another fre- 
quency order parameter. From Eq. H45I) . it is straight- 
forward to show D - {KAy^^ ^ Q^/^. One thus ex- 
pects D ~ iV-^/(2^) at criticality and D - Ar-i/4 
the detrained phase. We estimate from Fig. ^| that 
/3/(2P) = 0.105(10), which is also fully consistent with 
the estimate from Q. 
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FIG. 17: Critical decay of the frequency order parameter Q 
The straight line represents the best linear fit of slope 0.21. 
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FIG. 18: Critical decay of the frequency order parameter D. 
The straight line represents the best linear fit of slope 0.105. 



at finite values of K. 



VI. SUMMARY 

We have explored the collective behavior of locally cou- 
pled oscillators with random intrinsic frequencies on d- 
dimensional hypercubic lattices as well as globally cou- 
pled oscillators where each oscillator is coupled with ev- 
ery other one with the same strength. Both phase syn- 
chronization and frequency entrainment have been stud- 
ied. We have probed those phenomena through one typ- 
ical model of growing surface. By measuring the mean- 
square width for both the phase and the phase velocity 
of the growing interface, we have estimated some lin- 
earity characteristics of the system. In particular, the 
lower critical dimension for phase synchronization has 
been obtained to be four. For frequency entrainment, 
the lower critical dimension has been numerically found 
to be two. The effects of the nonlinear sine coupling 
on the phase synchronization as well as the frequency 
entrainment have also been investigated by means of nu- 
merical simulations; revealed is that the sine coupling 
tends to suppress not only phase synchronization but also 
frequency entrainment rather than enhancing those. It 
is found that phase synchronization emerges for d > 4 
and frequency entrainment transition occurs for d > 2. 
The nature of the phase synchronization transition in 
five dimensions has also been explored and the critical 
exponents (3 and i> have been measured. The values 
are observed to differ from those for the globally cou- 
pled system, manifesting that the five- dimensional sys- 
tem belongs to the different universality class from the 
mean-field one. It further allows the speculation that the 
upper critical dimension for phase synchronization might 
be higher than five, although further study is necessary 
for conclusive results. 



Finally, we comment on the possibility of the complete 
frequency entrainment (Q = 1) in the system of globally 
coupled oscillators. From Eq. (|46|l . one can easily notice 
that, a g{uj) has no cutoff at finite frequency u like our 
choice gitu) ^ exp(— a-'^/4cr), K should grow arbitrarily 
large to yield Q = 1. However, with finite cutoffs in g{u!), 
the system can exhibit complete frequency entrainment 



Acknowledgments 

One of us (MYC) was supported in part by the KOSEF 
Grant No. ROl-2002-000-00285-0 and by the BK21 Pro- 
gram. 



[1] For a list of references, see A.T. Winfree, The Geometry 
of Biological Time (Springer- Verlag, New York, 1980); J. 
Theor. Biol. 16, 15 (1967). 

[2] Y. Kuramoto, Chemical Oscillations, Waves, and Turbu- 
lence (Springer- Verlag, Berlin, 1984). 

[3] Y. Kuramoto, in Proceedings of the International Sympo- 
sium on Mathematical Problems in Theoretical Physics, 
edited by H. Araki (Springer- Verlag, New York, 1975); 
Y. Kuramoto and I. Nishikawa, J. Stat. Phys. 49, 569 
(1987). 



[4] G.B. Ermentrout and J. Rinzel, J. Math. Biol. 22 1 
(1985); H. Daido, Prog. Theor. Phys. 75, 1460 (1986); 
77, 622 (1987); J. Phys. A 20, L629 (1987); Phys. Rev. 
Lett. 68, 1073 (1992); H.-A. Tanaka, A.J. Lichtenberg, 
and S. Oishi, Phys. Rev. Lett. 78, 2104 (1997); H. Hong, 
M.Y. Choi, J. Yi, and K.-S. Soh, Phys. Rev. E 59, 353 
(1999); H. Hong, M.Y. Choi, B.-G. Yoon, K. Park, and 
K.-S. Soh, J. Phys. A 32, L9 (1999); J.A. Acebron, L.L. 
Bonilla, and R. Spigler, Phys. Rev. E 62, 3437 (2000); 
H. Hong and M.Y. Choi, ibid. 62, 6462 (2000). 



18 



[5] H. Sakaguchi, S. Shinomoto, and Y. Kuramoto, Prog. 
Thcor. Phys. 77, 1005 (1987); H. Sakaguchi, tbid. 79, 39 
(1988). 

[6] H. Daido, Phys. Rev. Lett. 61, 231 (1988). 
[7] M. Bahiana and M.S.C. Massunaga, Phys. Rev. E 49, 
R3558 (1994). 

[8] T. Aoyagi and Y. Kuramoto, Phys. Lett. A 155, 410 
(1991). 

[9] S.H. Strogatz and R.E. MiroUo, Physica D 31, 143 
(1988); J. Phys. A 21, L699 (1988). 

[10] See e.g., A. Pikovsky, M. Rosenblum, and J. Kurths, Syn- 
chronization: A Universal Concept in Nonlinear Science 
(Cambridge University Press, Cambridge, 2001). 

[11] H. Hong, Hyunggyu Park, and M.Y. Choi, Phys. Rev. E 
(in press). 

[12] S.F. Edwards and D.R. Wilkinson, Proc. Roy. See. Lon- 



don A 381, 17 (1982). 

[13] D. R. Nelson, Phys. Rev. B 26, 269 (1982); R. Cucrno 
and E. Moro, Phys. Rev. E 65, 016110 (2001). 

[14] S. Das Sarma and P. Tamborenea, Phys. Rev. Lett. 66, 
325 (1991); C. Herring, in The Physics of Powder Met- 
allurgy, edited by W. E. Kingston (McGraw-Hill, New 
York, 1951); W. W. MuUins, J. Appl. Phys. 28, 333 
(1957). 

[15] Sec, e.g., R.L. Burden and J.D. Fairos, Numerical Anal- 
ysis (Brooks/Cole, Pacific Grove, 1997), p. 280. 

[16] H. Hong, B. J. Kim, and M. Y. Choi, Phys. Rev. E 65, 
026139 (2002). 

[17] R. Botet, R. JuUien, and P. Pfeuty, Phys. Rev. Lett. 49, 
478 (1982). 



